
source("../src/communitySim.R")
source("../src/simCancer.R")

if(F) {
# how I got the 0.4899 figure
bob = function(qq) abs(qnorm(0.95, 0, sqrt(qq)) - log(10)/2)
optim(0.5, bob)
}
totalVar = 0.4899106
totalSD = sqrt(totalVar)
propGroupSD = 0.2



parameters = list(
	   geno=log(1.5), env=log(1.25), "geno:env"=log(2),
		 size=150000, probGeno=0.2,
		 probEnv=0.3, genoMiss= 0.1,
		 envMiss = 0.1, oneEnvPerCommunity =F,
     cancerMissRate=0.2, falseCancerRate=0.00045,
		 sdGroup=propGroupSD * totalSD,
		 sdPerson=sqrt(1-propGroupSD^2)*totalSD,
		 Enrollmentprobs = c(20000, 40000, 50000, 40000),
		 agerange = c(35, 69),
		 Followup = c(5, 10, 20, 30),
		 Ncommunity =  50,
		 EqualCommunity=T
)



sig = 0.001



populationData = getPopData() 
 
if (interest=="diabetes"){
deathFile=list(
          "F"=read.table("../data/MortalityFemale.txt",header=T),
          "M"=read.table("../data/MortalityMale.txt", header=T)
)


# common cancer means specific disease being selected. For cvd, we choose all of them  
CommonCancer=NULL


lambda =list(F=list(x=seq(30, 85, by=5), y= exp((-0.5)*0.585)*c(2.7, 3.9, 4.6, 6.6, 9.6, 12.3, 15.5, 18.2, 18.8, 
18.4, 16.4, 13.3)/1000),
M=list(x=seq(30,85, by=5), y=exp((-0.5)*0.585)*c(2.2, 3.7, 5.7, 8.5, 11.9, 15.8, 19.4, 23.3, 23.1, 22.3, 20.1,
14.7)/1000) )

DiffCancer=NULL
} 

if (interest="cancers"){

CommonCancer = c("Colon","Lung","Breast", "Prostate", "Stomach")

deathFile = list(
  "F"=read.table("../data/MortalityNocancerFemale.txt", header=T),
  "M"=read.table("../data/MortalityNocancerMale.txt", header=T)
)
DiffCancer= list(
	"F"=lambdaDiffCancer(myfile="../data/IncidentFemaleDiffCancer.csv",
	 	agegroup),
   M=lambdaDiffCancer(myfile="../data/IncidentMaleDiffCancer.csv",
	 	agegroup)
	 )
# change simcohort and getDeath functions to take deathData and DiffCancer arguments
# instead of reading files, so that it's easier to use cvd death data
# if DiffCancer is missing (it will be for CVD), make all event types the same.


lambda =list(F=list(x=seq(30, 90, by=5), y= exp((-0.5)*0.585)*c(116.4,169.7,272.6,401.8,555.4,746.8,1003.2,1254.8,1528.7,1741.4,
1903.1,1950.2,1922.5) / 100000)  ,
M=list(x=seq(30,90, by=5), y=exp((-0.5)*0.585)*c(61.8,81.2,123.6,247.8,466.8,877.5,1403.7,2096.6,2598.4,2887.2,
3086.4, 3156.3, 2927.5)/ 100000) )


}

if (interest="cvd"){
deathFile = list(
  "F"=read.table("../data/MortalityNoCardiFemale.txt", header=T), 
  "M"=read.table("../data/MortalityNoCardiMale.txt", header=T)
)

DiffCancer= list(
	"F"=lambdaDiffCVD(myfile="../data/IncidenceFemaleCVD.csv",
	 	agegroup),
   M=lambdaDiffCVD(myfile="../data/IncidenceMaleCVD.csv",
	 	agegroup)
	 )	
# common cancer means specific disease being selected. For cvd, we choose all of them 
# program will add "all" automatically 
CommonCancer=NULL

lambda =list(F=list(x=seq(35, 65, by=5), y= exp((-0.5)*0.585)*c(271.95, 343.05,463.46,658.36,798.68,1065.21, 1315.74) / 100000)  ,
M=list(x=seq(35,65, by=5), y=exp((-0.5)*0.585)*c(265.36,407.74,603.52, 846.10, 1077.14, 1407.05, 1649.21)/ 100000) )

}

Nsim=100



theGender=c("F","M")
Percentile = c(0.025, 0.975)

carray=array(NA, c(length(parameters$Followup),length(theGender),Nsim))
dimnames(carray)= list(as.character(parameters$Followup),theGender,NULL)


for(D in 1:Nsim) {
 carray[,,D] = CaseNum( parameters, lambda, CommonCancer,populationData, deathFile, theGender=c("F","M"))
                }

EstiCase= apply(carray, c(1,2), function(qq) mean(qq))

EstiCI = array(NA, c(length(parameters$Followup), length(theGender), length(Percentile)))
dimnames(EstiCI) = list(as.character(parameters$Followup),theGender, as.character(Percentile))

EstiCI = apply(carray, c(1,2), function(ci) (quantile(ci,probs=Percentile) ))

interNumCI=list(CI = EstiCI, cases = EstiCase)

FCI=round(cbind(as.matrix(interNumCI$cases[,1]),t(as.matrix(interNumCI$CI[,,1]))))
MCI=round(cbind(as.matrix(interNumCI$cases[,2]),t(as.matrix(interNumCI$CI[,,2]))))
interEstCI=cbind(FCI,MCI)
colnames(interEstCI)=c("F","2.5%","97.5%","M","2.5%","97.5%")


#####################################################################################################

CaseNum = function(parameters, lambda, CommonCancer, populationData, deathFile, theGender=c("F","M")){

  
  if(is.null(CommonCancer)){
CommonCancer=	c(colnames(DiffCancer[[1]]) [ !colnames(DiffCancer[[1]])== "Total" ],"All")
}   else {
CommonCancer = c(CommonCancer, "All")
}

datamat = simCohort(parameters, lambda, populationData, deathFile=deathFile, DiffCancer)

size=parameters$size
Followup=parameters$Followup

NumOfCases=array(NA, c(length(Followup),length(theGender)))
dimnames(NumOfCases)= list(as.character(Followup),theGender)

verbose=F
for(Dcancer in CommonCancer) {
	if(verbose) cat(Dcancer)
	thisCancer = rep(F, parameters$size)
# find the cancer we're interested in
  if (Dcancer=="All"){
  thisCancer = datamat$Cancer
  }else{                                  
	thisCancer[grep(Dcancer, datamat$CancerType)] = T
	}
	
	for(Dfollowup in Followup) {

	 	# set thisEventTime equal to age at end of followoup period
	 		ageEndFollowup = datamat$Age + Dfollowup - datamat$Enrollment
		hadEvent = datamat$Event < ageEndFollowup & thisCancer
	 	

		# eventType=T if the cancer of interest occurred during the followup period
		
    NumOfCases[as.character(Dfollowup),1]=sum(as.numeric(hadEvent&datamat$Gender=="F"))
    NumOfCases[as.character(Dfollowup),2]=sum(as.numeric(hadEvent&datamat$Gender=="M"))

		
   }

}

return(NumOfCases)
}

